Inhomogeneous spectral moment sum rules for the retarded 
Green function and self-energy of strongly correlated electrons or 
ultracold fermionic atoms in optical lattices 

J. K. FreerickJ*] 

Department of Physics, Georgetown University, Washington, D. C. 20057 

V. Turkowski 

Department of Physics and Nanoscience and Technology Center, 
University of Central Florida, Orlando, Florida 32816 
(Dated: July 26, 2009) 



1 



Abstract 

Spectral moment sum rules are presented for the inhomogeneous many-body problem described 
by the fermionic Falicov-Kimball or Hubbard models. These local sum rules allow for arbitrary 
hoppings, site energies, and interactions. They can be employed to quantify the accuracy of 
numerical solutions to the inhomogeneous many-body problem like strongly correlated multilayered 
devices, ultracold atoms in an optical lattice with a trap potential, strongly correlated systems 
that are disordered, or systems with nontrivial spatial ordering like a charge density wave or a spin 
density wave. We also show how the spectral moment sum rules determine the asymptotic behavior 
of the Green function, self-energy, and dynamical mean field, when applied to the dynamical mean- 
field theory solution of the many body problem. In particular, we illustrate in detail how one can 
dramatically reduce the number of Matsubara frequencies needed to solve the Falicov-Kimball 
model, while still retaining high precision, and we sketch how one can incorporate these results 
into Hirsch-Fye quantum Monte Carlo solvers for the Hubbard (or more complicated) models. 
Since the solution of inhomogeneous problems is significantly more time consuming than periodic 
systems, efficient use of these sum rules can provide a dramatic speed up in the computational 
time required to solve the many-body problem. We also discuss how these sum rules behave in 
nonequilibrium situations as well, where the Hamiltonian has explicit time dependence due to a 
driving field or due to the time-dependent change of a parameter like the interaction strength or 
the origin of the trap potential. 

PACS numbers: 71.27.+a, 71.10.Fd, 73.21.-b, 03.75.-b, 72.20.Ht 
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I. INTRODUCTION 



Spectral moments are integrals of powers of frequency multiplied by the corresponding 
spectral function and integrated over all frequency. They can reveal important information 
about the structure and spread of the spectral function and, in some cases, can also reveal 
interesting information about different many-body correlation functions. In theory, knowl- 
edge of all spectral moments allows one to reconstruct the function, but that procedure is 
well-known to be unstable and is not commonly used in numerical calculations. The spectral 
moments also correspond to derivatives of the Green functions with respect to relative time 
(either real time or imaginary time), evaluated at the point where the relative time is zero. 
As such, the spectral moments provide information about the relative-time dependence, 
when expanded as a Taylor series in time. 

Spectral moment sum rules for the many-body problem were investigated in 1967 by 
Harris and Lange^ shortly after the Hubbard model^ was introduced. In that work, one can 
find the moment sum rules for the first three moments of the retarded Green function of the 
Hubbard model and various approximations like the alloy analogy problem (which is equiv- 
alent to an inhomogeneous Falicov-Kimball model 3 ). They also developed a strong-coupling 
projection method to find spectral moments within each of the different Hubbard bands. 
This approach is an approximate one, as the projection operator is developed in a power se- 
ries of the hopping divided by the interaction strength. Shortly thereafter, Nolting 4 applied 
the spectral moment sum rules to develop approximations for the momentum-dependent 
Green function that have two poles, with the weights and locations of the poles fixed by 
the corresponding sum rules. This approach has been developed quite extensively, and ap- 
plied to a variety of different problems including dynamical mean-field theory^^^. That 
work examined ferromagnetic and antiferromagnetic long-range order in the Hubbard model, 
different approximation schemes for perturbation theory that produce the correct strong- 
coupling limit, and extended the sum rules to the retarded self-energy. The approach has 
also been applied to photoemission and inverse photoemission^, where it was recognized 
that the moments of the so-called lesser and greater Green functions play an important role. 
Finally, Harris and Lange's projection technique was applied to the lesser and greater Green 
functions to examine spectral properties of the Hubbard modet^. 

Most of that work has had as its focus using the sum rules to develop different approx- 
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imations, or to learn semiquantitative features of the many-body problem. But there is 
another application of spectral moment sum rules that is quite important for computational 
work. The sum rules allow, in certain circumstances, for numerically exact computations to 
be benchmarked against the sum rules. Steven White applied this to the Hubbard model 
in two dimensions^ 2 -, and Deisz, Hess, and Serene applied it to the fluctuation-exchange 
approximation^ at the Matsubara frequencies. The sum rules have been generalized to 
nonequilibrium situations for the Green functions 14 and self-energies^. In the case when 
one applies a spatially uniform, but time- dependent, electric field to the system, it turns 
out that many of the low-order moments are time-independent, even though the Hamilto- 
nian has explicit time- dependence due to the field being turned on at a specific time (or 
because the field has nontrivial time dependence). Nonequilibrium dynamical mean- field 
theoryi&i 7 -^ can solve the many-body problem exactly, and the sum rules are employed to 
benchmark the quality of the solutions 19 . 

Finally, spectral moment sum rules, when viewed as a Taylor series expansion in time, are 
now being employed to improve both the speed and the quality of Hirsch-Fye quantum Monte 
Carlo approaches 2 ^ for solving the impurity problem in dynamical mean-field theory 2 ^. In 
particular, recent work— has shown that employing the exact Taylor series expansions 
for short imaginary times allows one to use much larger Trotter error, yet achieve high 
accuracy with the Hirsch-Fye algorithm, so that it is competitive with continuous-time- 
based approaches 2 ^. 

In this work, we want to extend the spectral-moment sum rules for the retarded Green 
function and retarded self-energy to inhomogeneous cases, which are becoming more impor- 
tant, and where one can find similar improvement in the efficiency of impurity solvers for 
use in inhomogeneous dynamical mean-field theory problems. There are four main thrusts 
of work in the inhomogeneous many-body problem currently: (i) examining the properties 
of strongly correlated multilayers because of their potential for quantum-mechanical engi- 
neering of device properties 2 ^; (ii) examining the properties of ultracold atoms on optical 
lattices but spatially confined by a magnetic or optical trap 27 ; (iii) examining the proper- 
ties of a strongly correlated material that is disordered; and (iv) examining properties of 
strongly correlated materials with inhomogeneous spatial ordering (like a charge or spin 
density wave). In the multilayer problem, most solutions have relied on approximate tech- 
niques for the Hubbard model, or have examined simplified models like the Falicov-Kimball 
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model (recently, however, the numerical renormalization group has been applied to multi- 
layer Hubbard systems 2 ^ 2 ^) . In the cold atom problem, little work has been applied to Green 
function-based techniques (although this is increasing); here the spectral-moment sum rules 
could aid both in benchmarking and in improving the accuracy and efficiency of numerical 
algorithms. The strongly correlated material with disorder problem has had many different 
techniques applied to it, but most of them are approximate in one way or another, so un- 
derstanding the quality of the approximations is important. Less work has taken place in 
ordered phase calculations, especially for sum rules when the system is spatially ordered; 
these results can be immediately examined with the results presented here. In all cases, 
use of sum rules can help provide quantitative data to analyze how accurate the numerical 
solutions are. Finally, inhomogeneous systems are likely to be studied within the nonequilib- 
rium context. Here, we also generalize the nonequilibrium sum rules to the inhomogeneous 
environment and we consider a wide range of different nonequilibrium contexts. 

The remainder of this contribution is organized as follows: in Section II, we discuss the 
formalism for deriving the spectral moments in equilibrium; in Section III, we generalize to 
nonequilibrium cases; in Section IV, we discuss different applications of sum rules within 
dynamical mean-field theory; and in Section V, we present the conclusions and summary. 

II. FORMALISM FOR DERIVING SPECTRAL MOMENTS IN EQUILIBRIUM 

The equilibrium Hamiltonian for the inhomogeneous Falicov-Kimball 3 and Hubbard 2 
models can be written in the following unified form: 

n = X> c h- -Hiflfj + E MfA^ (i) 

ij ij i i i 

Here c\ (cj are the creation (annihilation) operators for a conduction electron at site % and 
f\ and f i are the corresponding operators for a localized electron at site i; in the case of the 
Hubbard model, the c-electrons are the spin-up electrons and the /-electrons are the spin- 
down electrons. The hopping matrix is denoted by —t^ and —t{- for the c- and /-electrons, 
respectively; for the Falicov-Kimball model we have v = 0, while for the Hubbard model 
we have v — t (one can also consider an asymmetric Hubbard model with < t* 7^ t). 
Note that the hopping matrix need not be translationally invariant, the only requirement 
is that it is Hermitian. The local chemical potentials are denoted by /Xj = /1 — Vi and 

5 



Hfi = [if — Vfi for the conduction and localized electrons with Vi and V/i the corresponding 
local potentials (/ij = for the Hubbard model in a vanishing magnetic field; /ifi — for 
the Falicov-Kimball model, since it does not enter into the conduction electron moments), 
and the local interaction between different particles is C/j. In the case of disorder, one often 
averages the Hamiltonian, and different measurable operator expectation values, over the 
given disorder distribution for the different parameters that are disordered. We will not 
discuss the details of how to treat those kinds of problems here. The sum rules we derive 
would be for one quenched disorder configuration (corresponding to a particular choice of 
parameters in the Hamiltonian), and could subsequently be averaged with respect to the 
given disorder distribution, if desired. Note that this Hamiltonian is time independent, and 
we have no net current flow, so it corresponds to an equilibrium problem (in other words, 
it is in the "slow limit" if the potentials correspond to an electric field, where the charge 
rearranges itself into a static redistribution in response to the potential rather than allowing 
current to flow). We will also examine a wide class of nonequilibrium cases below. 

We do not make any assumptions about the translational invariance of any of the pa- 
rameters that enter the Hamiltonian, so the Green function will generically depend on two 
spatial coordinates rather than on their difference. But in equilibrium, the system does have 
time-translation invariance, so we can describe the Green functions with a single frequency 
by making a temporal Fourier transform. The retarded Green function is defined to be 



in the time representation, with Z = Trexp(— (37i) the partition function, f3 the in- 
verse temperature, and {A, B} + = AB + BA the anticommutator. The creation and 
annihilation operators are in the (equilibrium) Heisenberg representation, where 0(t) = 
exp(iHt)0 exp(—iHt) for any operator 0\ in this representation, the time-translation in- 
variance is easy to show due to the invariance properties of the trace and the fact that the 
Hamiltonian commutes with itself. The frequency representation for the retarded Green 
function is 






(3) 



The spectral function is then defined to be 




hmG*(u), 



(4) 
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duu n A$(u) 



and the spectral moments become 

Similar to the homogeneous case^ 1 ^, one can easily show that 

/^(R,R,) = -21m 



r ^ G 3 (t ' 0) 



(5) 



(6) 



J t=o+ 



note that this representation uses the fact that 8(0) = 1/2, so the factor of 2 is needed in 
order to make 28(0) = 1 and have the correct normalization. Using the Heisenberg equation 
of motion, one can relate the time derivatives to commutators with the Hamiltonian, to 
obtain 

^(R,,R,) = Re({L"c,(0),c](0)} + ), (7) 

where L n O = [...[[0,Ti],Ti]...7{\ is the multiple commutation operator, performed n times. 
In Eq. ([7]), we used a shorthand notation, where the angular brackets denote the trace over 
all states weighted by the density matrix exp(— f3H,)/ Z. 

Evaluating the commutators is straightforward, but tedious. The results are 



/i^(Ri, Rj, T) — 5. 



/if (Rj, Rj, T) — —tij — 5ij(fii — UiUfi), 



(8) 
(9) 



//£(R», Rj, T) = ^ tukj + ilk - UiUfJUj + Ujdxj - UjUfj) 
1 

+ SiAbu - U,n ft ) 2 + Ufn fi (l - n fi )\, (10) 
/zf(Rj,Rj,T) = — / ] tjktkitij 

kl 

- (fjii - UjUf^yjutij - y]tu(jj,i - Uin f i)tij - ytuti^Hj - UjUfj) 

1 1 1 

- (ni - Uin fi )% - (ni - r,// /7 )/,;(//, - UjUfj) - tij(nj - UjUfj) 2 

- U 2 n fi (l - ///,)/,, - UiUjUj (fUifjfj) - ".fi".fj - tijUjn fj (l - n fj ) 

- Sij [(/Ji - UiTifif + 3Ufmn fi (l - n fi ) - Ufn fi (l - n fi )(l + n fi )] 

l,m 

+ %c/^(/if-//{) [tlifltt + tiiflfi)} +s ij u^t f u(flfi) 

I I 

- SijU^Ui [tliflficjci) +t{ l (fJf l 4c l )]+U i U j [/];(/;./-';/•,) • // ; (./y.//yv) 

1 

(11) 



where n/j = this result corrects a factor of 2 error in the third line of Eq. (29) of 



Ref . 1 1 5l. Note, that in cases where we have inhomogeneity arising from the spatial long-range 
order (say charge or spin density wave order, for example), then the Hamiltonian is actually 
translationally invariant, but the sum rules have inhomogeneous results due to the explicit 
evaluation of the different expectation values (the charge density will vary from site to site 
in a charge density wave, for example). This then gives rise to spatially inhomogeneous 
moments. 

Next, we examine the retarded self-energy moments. To begin, we start with the Dyson 
equation in the real space for the frequency- dependent Green function and self-energy: 

GfM = GfV) +J2 G TH^)G^), (12) 

ki 

where GfP(u) is the noninteracting retarded Green function on the lattice. When the 
frequency uj is large enough, all Green functions and self-energies become purely real (on 
the infinite-dimensional hypercubic lattice, there can be an exponentially small imaginary 
part since the bandwidth is infinite, but this plays little role for large frequencies), and by 
using the spectral formulas, 

G « M ._ir^, (13) 

J 7T J_ 00 UJ - Uj' + lb 

SgH = --/ dJ ^- + £^ = oc), (14) 

7T J _ 00 UJ — Uj' + ZO J 

we can expand the functional dependence of the Green function and self-energies in terms 
of the corresponding moments (since the denominators in the integrals never vanish when 
the numerators are nonzero) yielding 



m=0 



E,» = E,> = oo) + £^M, (16) 

m=0 



where 



1 f°° 

C^R h Rj) = — / dujuj^m^uj), (17) 

are the spectral moments for the self-energy [S^(u; = oo) is a real constant equal to the 
large-frequency limit of the self-energy] . The expansions in Eqs. ([15]) and ffl6|) (and a similar 



-mjR-j)) (20) 
mi Rj) 



expansion for the noninteracting retarded Green function) are substituted into the Dyson 
equation in Eq. f|T2|) and then one equates powers of 1/ui to find 

uf(Rj,Rj) = /tf(Rj,Rj) (18) 
/if (Ri, Rj) = /if (R, Rj) + /tf(Ri, R)SL(^ = oo)/if (R m , Rj), (19) 
/if (Rj, Rj) = /if (Rj, Rj) + /if (Rj, R i )Sf l (u; = oo)/if (R m , Rj) 

+/if (Rj, Ri)C(f (Rj, R m )/if (R m , Rj) 

+/if(Rj,R)£f> = oo)/if(R m ,R,), 
/if (Rj, Rj) = /if (Rj, Rj) + /if (Rj, R)Sfj(^ = oo)/if (R, 

+/if (Rj, R;)Cf (R;, R m )/if (R m , Rj) 

+/if (Rj, R;)Cf (Rj, R m )/if (R m , Rj) 

+/if(R,R)Sf m (^ = oo)/if(R m ,Rj) 

+/if (Rj, Hi)Cq (Rj, R m )/if (R m , Rj) 

+/if (R, R*)E£> = oo)/if (R m , Rj), (21) 

where the matrix /if (Rj, Rj) is the nth spectral moment of the noninteracting retarded 
Green function on the lattice. Those noninteracting moments are found from Eqs. (TTTj) 
with Ui = 0. Substituting the explicit values of the moments into Eqs. (fT8]) - (l2T]) finally 
yields the self-energy moments: 

Eg (u = oo) = SyUiUfi, (22) 

C^R,^) = dijUfrifiil - n fi ). (23) 

Note that the algebra required to arrive at these results is nontrivial. In cases where i and j 
are farther apart than the range of the hopping matrix (which is often chosen to be nonzero 
only for nearest neighbors) many moments are identically zero, but nontrivial cancellations 
are required to ensure that all of the off-diagonal moments vanish (no local approximation 
has been made for the self-energy here — these results hold in all dimensions). The first 
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self-energy moment is: 



Cf (R, Rj) = %r,-// /; (l - n,*) [£/<(! - n/i) - 



l,m 

i i 



+UiU j \1^ i (f}f i clc i )+t( j {f}f j clc i ) 



(24) 



which contains an off-diagonal term when one is in finite dimensions (it is local in infinite 
dimensions due to the scaling of the hopping matrix element with dimension) . 

The expression for the nonhomogeneous Green function moments Eqs. (l8l)- ffTTl) and for 
the self-energy moments Eqs. fl22l) - fl24l) are rather general. In particular they can be used to 
evaluate the moments in the case of a particular (quenched) configuration of disorder. Then 
we would average over some disorder distribution [say, P(Vi) for diagonal disorder]. This 
procedure requires one to perform calculations over a range of different disorder distributions 
and then perform the averaging. If the system tends to self-average, not too many specific 
configurations would be needed, but if there is interesting physics arising from rare regions 
of the distributions, many calculations would be needed, and these calculations can get to 
be rather lengthy. 



III. GENERALIZATION TO NONEQUILIBRIUM SITUATIONS 

One rather general form for the nonequilibrium Hubbard- Falicov-Kimball Hamiltonian is 

m = ■■^i^iy,r l ■■Y.ii/n.rii) ^/M/H- >>,;(/)/,./; 

ij ij i i 

i 

In this case, we have added time dependence to all of the parameters in the Hamiltonian, 
but have not introduced any additional forms of interaction within the Hamiltonian. Nev- 
ertheless, this generalization allows for a rather rich class of nonequilibrium problems to 
be studied. For example, if we are examining a multilayered device with electronic charge 
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reconstruction 3 ^ (where the potentials Vi and V/ are determined by an additional semiclas- 
sical Poisson equation), and we use a vector potential to describe an external electric field 
that drives current through the system 2 ^^, then we would have a time dependent hopping 
determined by the Peierls substitution 32 . If we want to examine an interaction quench, as is 
often studied in cold atom systems, we would have a (typically harmonic) trapping potential 
Vi and V/ and the interaction U% would become time dependent switching from one value for 
early times to another value for later times^ 3 - such as would occur near a Feshbach resonance 
if the bias magnetic field is changed from one value to another (for some experiments, the 
potentials or hopping could also change when the coupling changes); the switching could be 
sudden as in a rapid quench, or adiabatic, with a slowly varying change, or anything in be- 
tween. In addition, within the cold-atom picture, we could imagine creating time-dependent 
trap potentials Vi(t) and V/(t). This would allow us to examine what would happen if we 
applied an impulse to the atomic cloud, or if we shifted the origin of the harmonic potential 
from one spatial location to another, and then examined how the center-of-mass oscillates 
and damps back to the thermal state, or to some nonthermal steady state. Finally, we could 
examine the so-called Bragg spectroscopy experiment, where the optical lattice potential 
amplitude is oscillated with some set frequency and one observes things like the change in 
the momentum distribution after the system is driven for a certain period of time, or the 
change in the double occupancy. In this case, the hopping, the local potentials, and the 
interaction could all become time dependent, but can still be described by the general form 
of our Hamiltonian. 

There are a few subtle points to keep in mind with these nonequilibrium problems. For 
example, in the case of a multilayered system with electronic charge reconstruction, the 
charge reconstruction is created in the distant past, so it corresponds to an equilibrium static 
potential which does not contribute to any flow of current. The field that drives the current 
is described by a time-dependent vector potential, which is turned on at a particular time. 
One can examine the transient current flow or the steady-state current flow by examining 
how the system responds to the external field. 

Given this general form for the Hamiltonian, we next need to derive the sum rules. 
We work in a Heisenberg picture, because the operator algebra for the time-dependent 
creation and annihilation operators, at equal times, is unchanged from the standard fermionic 
commutation relations. The only difference from the nonequilibrium derivations worked out 
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previously^ 1 ^ is that here we need to work in real space for all calculations, because there 
is no translational invariance. In order to evaluate the expression for the nonequilibrium 
spectral moments, it is convenient to introduce the relative t = t\ — and the average 
T — (ti + t 2 )/2 time coordinates for the retarded Green function Eq. (|2J). In this case, 
the physical time at which one would likes to calculate the moments will correspond to 
the average time T, and the Fourier transform to frequency space must be performed with 
respect to the relative time t. Then, one can define the nonequilibrium Green function 
moments by generalizing the expression in Eq. (jSJ): 

1 f°° 

fM^(Ri, Rj, T) = / duu n ImG*(T,uj). (26) 

" J — oo 

In a similar way, one can define the nonequilibrium moments for the self-energy: 

C^(R 4 ,R„T) = — / duu m hnE*(T,u) (27) 

(for more details, see Ref. 15). By using these equations, one can easily show that the 
expressions for the nonequilibrium retarded Green function moments in Eqs. (|8l)-(fTTl) re- 
main unchanged in the nonequilibrium case, except the model parameters and the operator 
expectation values are replaced by their time-dependent forms: Ui — > L^(T), ^ — > yUj(T) 
and nfi — > rifi(T). Also, the operators in the correlation functions in Eq. ffTTl) are in the 
Heisenberg representation at the average time T, at which the spectral moment is calculated. 
Similarly, one can show that the expressions for the large frequency self-energy in Eq. ( 122]) 
and for the zeroth self-energy moment in Eq. ( 1231) remain the same in the nonequilibrium 
case. However, the expression for the first moment in Eq. ( 1241) will acquire an additional 
term proportional to the second derivative of the large frequency limit of the self-energy 
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Y,ij(T,u = oo). Namely, in the nonequilibrium case one finds: 
CftR^T) = 5ijU^(T)nfi(T)[l — n/j(T)][[/j(T){l — n/i(T)} — /i;.(T)] 



+<WCO£ [&(^*LCO<tf(mCO> +t{ m (T)tUT)(f i (T)f l (T)) 

l,m 

-2t{ m {T%{T)(f}{T)f m {T)) 

+8mt) E (W ( T ) - rfcn) kco<//co/iCo> + tkT){fhT)h{T)) 



i 

+ 4(^)(//(^)/KT)cJ(T) Q (T))^ 
UWU^UjiTWl^U^fjiT)^)) - n fi (T)n fj (T)] 
+U i (T)U j (T) [t f ^T)(f}{T)ttT)c){T)c % {T)) 

+ t{AT)(fKT)f J (T)c](T)c l (T)} . 



(28 



For most nonequilibrium problems we would consider, except for an adiabatically changing 
interaction, or an amplitude oscillation of the optical lattice, the derivative term would 
vanish almost everywhere. 



IV. APPLICATION OF SPECTRAL MOMENT SUM RULES TO DYNAMICAL 
MEAN-FIELD THEORY 

There are two immediate applications of spectral moment sum rules within dynamical 
mean-field theory. The first is to use them to evaluate the high-frequency asymptotic be- 
havior exactly and then supplement the high frequency results by numerical calculations 
at low frequencies, and the second, closely related, is to use them to evaluate the short 
imaginary time behavior exactly and supplement with long-time numerical calculations. 
The latter has been already discussed within the context of the Hirsch-Fye quantum Monte 
Carlo algorithm for solving the impurity problem in dynamical mean-field theory and the 
results there show great promise as a means to improve the accuracy and the efficiency of 
calculations^ 1 ^. The basic idea (at half-filling) is that the curvature, which grows sharply 
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with increasing U, immediately determines the short-time behavior for the Green function. 
Because of the sharp dependence on r, one would need to use a very small discretization 
step for the QMC to accurately describe such behavior, which would then be very costly in 
terms of computational time. Instead, one uses a coarser grid, but before performing the 
Fourier transform to Matsubara frequencies, one simply creates a finer grid and uses the 
short time relations to find the behavior close to r = 0, and uses simple interpolation for 
other time values (a shape preserving spline would work well in this context). Then the 
Fourier transformation will much more accurately reflect the true behavior of the system, 
and all of the high-frequency structure will be properly recovered, so that the QMC can 
be used to determine the low-frequency data where it is most accurate. (It has already 
been demonstrated that this approach is competitive with other QMC techniques such as 
the continuous-time algorithm.) Issues of accuracy and efficiency will become increasingly 
important for inhomogeneous dynamical mean-field theory problems (like multilayers or ul- 
tracold atoms in a trap) because one needs to solve an impurity problem at each inequivalent 
lattice site of the inhomogeneous system. This can range from tens to hundreds of impurity 
solvers for multilayered systems to many thousands or more for ultracold atomic systems in 
a trap. We won't discuss the application within quantum Monte Carlo approaches further 
here, and instead will concentrate on examining a different application, which is to solve the 
dynamical mean-field theory for the Falicov-Kimball model with fixed local chemical poten- 
tials on the lattice sites. This approach works equally well (with appropriate modifications) 
for the two-site approximation to the Hubbard model^ in the insulating phase. 

The dynamical mean-field theory for the Falicov-Kimball model is well established in the 
literature^ 1 ^, so we just present the relevant formulas. Starting from a local self-energy 
Ej(io; n ), one must solve the Dyson equation for the local Green function [Gu(iui n )] 



Here u n = Tr(2n + l)//3 is the fermionic Matsubara frequency. When the system is homo- 
geneous, a Fourier transformation allows this problem to be solved immediately. When one 
has inhomogeneity in one dimension, as in multilayered structures, the local Green function 
(of an infinite device) can be found from the so-called quantum zipper algorithrn^L^ which 
expresses the Green function in terms of continued fractions. For finite cold atom systems, 
one can use LAPACK (or sparse matrix) routines to perform numerical matrix inversions 




(29) 



k 



14 



to find the local Green function^^. Whatever the technique, we assume that one can solve 
the Dyson equation to determine the local Green function. Next, the effective medium 
and dynamical mean field \ are extracted via the scalar equations 

G^lUn) = {[GtiiiLOn)}- 1 + T.iilUn)}' 1 . (30) 

and 

Xi{iu n ) = iu n +Hi- [Gifan)]' 1 - (31) 

[All quantities in Eqs. fl30l and fl3Tl) are scalar quantities — there are no matrix operations 
here.] Then one calculates the filling of the local electrons nfi = Zu/(Zoi + Zu) with 



Z Ql = e^' 2 TT (32 ) 
- LJ - iu n 

71= — CO 

and 



n=—oo 

Now we find the new local Green function 

Gn(iu} n ) = - — \ + - — tr~\ — FT' ( 34 ) 

and finally extract the local self-energy 

£*(«*>«) = [G-i^nT 1 - G« x («*>«)• (35) 

These equations are then iterated until they converge. The conduction electron filling sat- 
isfies 

Pd = -j ^2 G ™( iuJ n ); (36) 

n=— oo 

this result is not part of the dynamical mean-field theory iteration, but it is needed if one 
wants to update the chemical potential during the iterations to achieve a particular electron 
filling. Note that this summation is ill-defined and needs to be properly regularized (see 
below) . 

Typically, one chooses a set number of Matsubara frequencies, usually with an energy 
cutoff many multiples of the noninteracting electron bandwidth, and solves the (now finite) 
set of equations for the Green functions and self-energies by iteration starting from E, = 0. 
One can try to include the effects of the neglected tails of the summations and infinite 

15 



products to improve the accuracy and minimize the effect of the energy cutoff. By employing 
the exact sum rules, one can make this procedure work well. We describe this process next. 

First recall that we have already proven that the Matsubara frequency Green function 
and self-energy satisfy 

/4*j(Ri, Ri 

m=0 

and 



= £^3*. (37) 



M**n) = Si(oo) + £ (38) 

m=0 ^ n > 

when the Matsubara frequency is large \ui n \ ^> \U ma _^\+W in t, where W int is the half bandwidth 
(in real frequency) of the interacting density of states (valid only for finite-dimensional 
systems). These results follow from the definition of the spectral moments, and the spectral 
formula for the Green function and self-energy. Using these two relations, substituting into 
Eq. (13T)|) . and recalling the definition in Eq. f[3"Tj) . produces 

H^n) = -x- — + n—T- — \T~ + ° \T- — \3 ' 39 

2iu n 2 (iuj n y \{iu n yj 

This asymptotic expansion along with the expansion in Eq. f[3T|) will allow us to treat the 
tails in the summation for the conduction-electron filling in Eq. ( 1361) and in the infinite 
products needed for the localized-electron filling in Eqs. (1321) and (|33|) . 

We imagine taking an energy cutoff E c which is larger than the interacting density of 
states half bandwidth. Using that cutoff to determine the explicit Matsubara frequencies 
solved in the numerical implementation of dynamical mean-field theory, we examine only 
Matsubara frequencies with \u n \ < E c . This defines a cutoff integer n c corresponding to the 
Matsubara frequency closest to E c but lying below it. The tail for the conduction electron 
filling 5^I^ _2 Gu{iuj n ) / (3 + X^+i Gn(iuj n ) / (3 can now be evaluated analytically for the first 
four moments that we have calculated. Define the approximate Green function via 

riapprox/- \ 1 . ^f(Rt) R-i) /^^(Ri, Rj) /^^(Rj) Ri) ( ac\\ 

(j-jj \}W n ) — ~. 1 t- v5 1 T- v5 1 T- vZ — ■ l 4U / 

iu n (tu n ) 2 (tuj n y (%uj n y 

Then, the identity 

1 1 °° 1 

1 + exp (-Pfi) (3 iu n + fi 

allows us to evaluate infinite sums of inverse powers of iu n noting that the sum of l/iu n 
requires special regularization which is given by the result in Eq. fj4T|) . By taking derivatives, 
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and evaluating at \x = 0, we immediately learn that 

= /(») = 5. 



lUJr 



-y— 

if 1 



-/'(o) = - 







-J" (0) = 0, 

1 /3 3 
6 w 48 



Substituting into the expression for the conduction-electron filling then yields 

1/9 /? 3 1 nc 



48' 



yapprox 



(42) 
(43) 
(44) 
(45) 



(iw n )]. (46) 



n=—n c —l 

The summation, which needs to be evaluated numerically, vanishes rapidly for large o; n , so 
the filling can be computed quite accurately. 

Now we show how to evaluate the infinite products by properly taking into account the 
asymptotic limits using the spectral sum rules. Substituting the results from Eq. (J39l) into 
Eqs. ( 132]) and ( 1331) allow for the local electron filling to be computed. First note that 



n=n c +l 



Z 0l = e^ J] 

n 

and 

Zl . = e /3(/^)/2 e ^/ 



IH _ 1 1 _ 1 Uin fi - fij 
iuj ri 2(iu n ) 2 2 



2 n c 

n 

n=0 



1 



g?(iu} n )iu; r , 



(47) 



n 



n=n c +l 



1 + 



Hi — Ui 1 1 



2W, 



1 ?7tW/t - ^ 

2 (^ n ) 2 2 



2 n c 

n 

n=0 



ft 



1 

(48) 



The infinite products that have an infinite number of terms are approximated by rewriting 
the infinite product as the exponential of the sum of the logarithm of the individual terms. 
Replacing the sum by an integral (valid when the temperature is much smaller than the 
interacting half bandwidth) and converting the integral over frequency to an integral over 
z — 1/uj yields 

7 2 



3 f 1/Ec dz 

f / ^ln({l + te 2 } 2 + ^{a-c^} 2 ) 



n 1 + — + -. no + - — — ~ ° x p 

n=n c +l 

9) 

where we use a = \ii for Z 0i , a = fii — Ui for Zu, b = —1/2 for both, and c = — (L^n/, — /ij)/2 
for both. This then allows for an accurate evaluation of the filling using just the small set 
of numerical data generated for \n\ < n c . 
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FIG. 1: False grayscale image of the relative error of the (left) zeroth self-energy moment and the 
(right) third spectral function moment for a trapped atoms system confined in harmonic traps with 
a length scale of 30 lattice units, U = 5, = 0.2 and 51 x 51 lattice sites. Note that the scale 
for the figures is less than the full range of deviations to pick up the fine structure far from the 
center of the lattice. At the center of the lattice, the error is maximal and equal to approximately 
1.5% for the zeroth self-energy moment and approximately 0.01% for the spectral function. This 
is elaborated on in the next figure. 

Note that our use of the asymptotic expressions for the different many-body functions 
aided us in reducing the effort of computation only for problems that can be solved along 
the imaginary axis. This includes the Falicov-Kimball model (for static properties) and 
Hirsch-Fye quantum Monte Carlo techniques for other models (like the Hubbard model). 
There does not appear to be any simple use of these sum rules within real-frequency-based 
approaches like the numerical renormalization group. Of course, all of these results can also 
be applied to the homogeneous case. 

In Fig. [U we plot false grayscale images of two moment sum rules for each lattice site of 
a 51 x 51 square lattice. This is a system close to phase separation with parameters U = 5, 
T = 0.2, and a harmonic trap with a characteristic length scale of 30 lattice spaces for both 
the light and the heavy particles. Note how the errors, on the whole, are rather small. In the 
imaginary- axis calculation, which is used to determine the chemical potentials so that we 
have approximately 625 light and 625 heavy atoms on the lattice, we employ the use of the 
moment sum rules to sum the tails of the Matsubara frequency calculations, which reduces 
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FIG. 2: Density of states (a) and — Im£j(u;)/7r (b) for the central site of the lattice with the 
same parameters as the preceeding figure. Note how the DOS is smooth and rather featureless, so 
the sum rules work to high accuracy, while the self-energy has sharp features, which lead to less 
accurate moments. In the inset of panel (b), we show the large peak in the self-energy; the data 
with the solid line has a small step size of 0.0003i, while the dashed lines and the circles correspond 
to a step size of 0.004t. Note how the smaller step size is smoother. 

the computational time by about a factor of 10 for the same accuracy. In this real axis 
calculation, we use a frequency grid with a step size of 0.004t that runs from — 9.Qt to 9.61 
This step size does not allow us to pick up fine structure smaller than the step size. It also 
cannot pick up spectral weight lying outside of the bounds. Both of these issues can cause 
inaccuracies in calculations, especially when a system begins to order or phase separate; 
they don't enter too significantly for this example though; when we repeat the calculations 
with a smaller step size of 0.0003t, we find the error in the zeroth self-energy moment is 
reduced from 1.5% to 0.025% for the central site of the lattice. The temperature we use 
here is high enough that there is no order, and the calculations are under good control. In 
any case, we show the spectral function and — Im£j(a>)/7r for the central site of the lattice 
in Fig. [21 where the relative error for the zeroth moment of the self-energy is about 1.5%, 
and the relative error for the Green function is about 0.01%. Notice the sharp peak in the 
self-energy which leads to the higher error, while the Green function is quite smooth, and 
hence has very small errors. As the step size is reduced, the errors in the moments are also 
reduced, indicating that these errors are arising predominantly from the discretization size 
of the real frequency axis and the structure of the sharp features in the functions. 
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FIG. 3: (Color online) plot of the absolute error of the (a) zeroth self-energy moment and the 
(b) third spectral function moment for a multilayered inhomogeneous system described by the 
Falicov-Kimball model with electronic charge reconstruction. There are 30 self-consistent ballistic 
metal planes to the right and the left of the 20 plane thick barrier which has U = 6 and half-filling 
for the heavy electrons. The temperature is 1/(3 = 0.1, and the shift in the band centers AEf 
labels the different figures. Note how the errors for the self-energy are larger, while for the Green 
function the errors are small. This is because our grid size is too coarse to properly pick up the 
weight of the narrow peak in the self-energy. The relative error for the third moment of the DOS 
is less than 0.1% in the barrier; in the metallic leads, the third moment gets very small, and the 
absolute error arises primarily from the discretization of the numerical quadrature. The dashed 
line indicates the interface between the metal and the insulating barrier. 

We also examine the sum rules for inhomogeneous multilayered systems. In cases where 
there is no electronic charge reconstruction, we have found that our data satisfies all of the 
sum rules to high accuracy (typically better than 0.01%) except for cases with an insulating 
barrier where the self-energy develops a sharp peak at low frequency. Our frequency grid in 
previous calculations was sometimes too large to properly extract the zeroth moment sum 
rule for the self-energy, and errors could become very large because the numerical quadrature 
is greatly overestimating the weight within the sharp peak near u — 0. In cases where the 
peak is not so sharp, we once again find excellent agreement. A more challenging case, 
though, is a case where there is an electronic charge reconstruction, because the calculations 
become much more difficult numerically in this case, and we usually need to introduce a 
finite broadening into the calculation to be able to estimate the local DOS on each plane. 
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Hence, it is much more interesting to examine these cases for calculations of the sum rules. 
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FIG. 4: Density of states (a) and — Im£j(u;)/7r (b) for the plane with the largest error (plane 
number 37 with AE'j = 1) in the preceeding figure. Note how the DOS is once again smooth and 
rather featureless, so the sum rules work to high accuracy, while the self-energy has a sharp peak 
(maximum amplitude is 8000), which leads to less accurate moments. 

In Fig. [31 we plot the absolute errors of a self-energy moment and a spectral moment 
(the other self-energy moment appears similar, while the other spectral moments had much 
smaller errors). The system consists of a semi-infinite bulk ballistic metal attached to a 
sandwich of 30 ballistic metal planes, 20 Falicov-Kimball model planes and 30 ballistic 
metal planes, so the calculations are always for a thermodynamic limit system. Both the 
metallic leads and the barrier are at half filling, with a common chemical potential. We 
shift the center of the band of the barrier by the amount AEf and solve for the electronic 
charge reconstruction with a screening length of a few lattice spacings and a temperature 
of 1/(3 = 0.1. The imaginary axis solver did not use the summation of the tails, as we are 



using old data from Ref. |31|. The real-axis solver worked with a grid step size of O.Olt and 
ranged uj from — lit to lit. One can see that, similar to the cold atom example above, here 
we also see errors which are much larger for the self-energy than for the spectral function 
(we only show the barrier planes for the self-energy moment, since the self-energy vanishes 
in the metal). This arises for the same reason as before, but is more acute here, since the 
step size in frequency is larger. This is shown in Fig. H] for the plane with the largest error 
in the self-energy. Note how once again the DOS is smooth, which is why the moments are 
so accurate, but the self-energy has a narrow peak, whose weight is overestimated with the 
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coarse grid used in the calculation. 




FIG. 5: Density of states for the A sublattice at U = 1.5 with an electric field of strength E = 1 
turned on at time t = 0. The different panels correspond to different average times. Note how the 
main structure of the DOS in equilibrium, which consists of the singular peak and the finite peak 
is modified at the odd Bloch frequencies here to create additional structures that look reminiscent 
of the DOS of an ordered system (one can see small peaks near u = ±3 too). Modifications at the 
even Bloch frequencies can only be seen at short times. 

Spectral moment sum rules for ordered phases, like a charge-density-wave phase, have 
already been performed in equilibrium^ for the Falicov-Kimball model, and agree with the 
exact results to high accuracy. The spectral moment sum rules for the nonequilibrium 
case in a homogeneous system with a uniform electric field have also been verified to high 
accuracy^ 1 ^. In addition, in the rapid quench work of Ref. |33|, the retarded Green function 
is given by the equilibrium Green function of the particular value of the interaction for 
each average time. Hence, the moments, which hold in equilibrium, continue to hold in 
nonequilibrium with an interaction quench. Here we examine a nonequilibrium case at 

22 



half filling for both localized and itinerant electrons with charge density wave order in the 
presence of a uniform electric field at zero temperature. This system can also be solved 
exactly within dynamical mean-field theory, because the self-energy vanishes or is equal 
to U and it has no damping (i. e., it is real when expressed in the frequency-average time 
representation). Details of that work will appear elsewhere^ and follow closely the derivation 
of the nonequilibrium Green function for noninteracting electrons on a lattice 4 ^, but with 
some added complications due to the need for time-ordered products because of the CDW 
order. The local retarded Green function at half-filling with \x = U/2 satisfies 41 



G*(ti, h) = -i0(ti - t 2 ) ex P [^n( k > *i, h) + W 22 (k, t x ,t 2 ) ± Wi 2 (k, ti, t 2 ) ± W 2 i(k, h, t 2 )\ , 



where the sum over momentum is over the CDW Brillouin zone, which satisfies < 0, and 
the plus sign is for i e A sublattice and the minus sign for i e B sublattice. The time 
evolution operator U is a time-ordered product 



band structure and A(t) the vector potential, which is turned on at time t — [A(£) = 
— 0(t)Et]. The electric field is chosen to lie along the diagonal direction. The time-ordered 
product can be calculated directly for numerical work by employing the Trotter formula on 
the corresponding 2x2 matrices. 

Now we report on the nonequilibrium sum rules for the CDW phase. In this case the 
moments of the local Green functions on each sublattice satisfy the following: (i) the zeroth 
moment is 1; (2) the first moment is ±U/2; (iii) the second moment is 1/2+U 2 /4; and (iv) the 
third moment is ±[//4±[/ 3 /8 (the plus or minus signs correspond to the A or B sublattice). 
While, in theory, one can numerically calculate these moments by first finding the Green 
functions as functions of time, converting to average and relative time, Fourier transforming 
the relative time to a frequency, and finally calculating the moment sum rule by integrating 
over frequency (see below), this approach runs into a number of serious challenges. The most 
critical one is that the equilibrium DOS, when there is no field present, has a divergence in 
it that goes like an inverse square root of frequency. Such a singularity requires an infinite 
time domain to properly find the Fourier transform (because the function has an amplitude 



k 



(50) 




(51) 



We used the Peierls' substituted band structure with e k = lim^oot* Ylt=i cos(kja)/\/rf the 
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FIG. 6: Density of states for the A sublattice at U = 3 with an electric field of strength E = 1 
turned on at time t = 0. The different panels correspond to different average times. Note how the 
main structure of the DOS in equilibrium, which consists of the singular peak and the finite peak 
is modified at the even Bloch frequencies here to create additional structures that look reminiscent 
of the DOS of an ordered system (small peaks are near uj = too). Modifications at the odd Bloch 
frequencies can only be seen at short times. 

that decays like a power law in time), but our numerical calculations are always truncated to 
a finite range, so the Fourier transform has the singularity smoothed over and the truncation 
can lead to unphysical oscillations in the DOS as a function of frequency (due to the presence 
of a sharp cutoff in time). Furthermore, we calculate on a discrete grid in time, which can 
have further effects on the Fourier transform, especially for high enough frequencies. Finally, 
the results, particularly at large times, are sensitive to the number of energy points in the 
integration grid over the two energies. All of these challenges make it much more useful 
to directly calculate the results for the moments by evaluating the derivatives of the Green 
functions as functions of time. This can actually be done analytically for the form of the 
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Green function (because of the time-ordered products), and it produces exactly the requisite 
moments. In addition, we can do it numerically in the time domain, and we verify that the 
sum rules are satisfied to an accuracy much smaller than the step size in time (when evaluated 
in the time domain via numerical differentiation). 

Even though the numerical calculation of the DOS in the CDW phase is challenging, as 
explained above, we present results for this calculation in Figs. [5] and [6] for the A sublattice 
DOS. The parameters chosen are an electric field equal to 1, and turned on at time t — 0, 
and interaction strength U = 1.5 and U = 3 (we have a time domain cutoff running from 

— 500 to 500 with a time step of 0.0025). We compare a number of different average time 
results including (a) the equilibrium result t ave . — > — oo, (b) t ave = 0, (c) t ave = 10, and (d) 
^ave. = 500. We do not use the "equilibrium" results from the real time calculation because 
our finite time-domain cutoff leads to spurious oscillations. In addition, the DOS is already 
modified at t av0 . = because the relative time Fourier transform involves a large number 
of points in time where the field is on, and because the Green function has structure with 
such long-ranged tails in time, one can see an effect even before the average time where the 
field has been turned on. When the field is turned on, the DOS naturally changes shape 
(and the square- root singularity appears to be smoothed into a finite peak), but the changes 
are much smaller than in the normal phase. Surprisingly, we predominantly see structure 
at either odd Bloch frequencies or even Bloch frequencies, but not both, and the structure 
does not look like a broadened, and split delta function as in the normal phase, but instead 
looks more like some kind of ordered phase gap structure, but the features can be quite 
small in some cases. It is clear we have not yet fully reached the steady state, but it also 
appears clear what the steady state DOS will eventually look like. Finally, we comment that 
if we take the DOS as functions of frequency (our frequency range is chosen to run from 

— 10 < uj < 10), multiply by the appropriate power of frequency and integrate to check the 
sum rules, then all sum rules except the second moment are satisfied to better than 1%. 
The second moment is worse, because we get significant contributions from the noisy tails 
of the DOS which don't cancel as they do for the odd moments. If we put the frequency 
cutoff closer to \u\ < 5, then the second moment accuracy increases dramatically. But we 
know that this is the least accurate way to check the moment sum rules. 
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V. CONCLUSIONS 



In this work, we have shown how to determine exact spectral moment sum rules for the 
retarded Green functions and self-energies of inhomogeneous strongly correlated systems. 
While our analysis is quite general, for concreteness, we provided explicit results for a com- 
bined Hubbard-Falicov-Kimball model. We envision these results can be applied to strongly 
correlated multilayers, ultracold atomic systems in traps, or strongly correlated materials 
with disorder (although we did not discuss the last case in much detail). The sum rules can 
be used to gain qualitative information about the many-body solutions, benchmark numeri- 
cal calculations, be used to evaluate the tails of infinite sums (or infinite products) allowing 
for smaller energy cutoffs, or be used to improve the accuracy of Hirsch-Fye quantum Monte 
Carlo approaches by determining the short-imaginary-time behavior exactly. We provided 
a full derivation of the sum rules in equilibrium, and then discussed a number of different 
nonequilibrium situations appropriate for multilayers, for cold atom systems, and for ordered 
phase systems (such as a CDW). In the case of multilayers, the change in the local electrical 
potential energy, induced by an electronic charge reconstruction, modifies the sum rules, 
but the additional scalar potential, that creates an electric field to drive current through the 
device, does not. This motivates one to decouple the description, using a scalar potential 
and Poisson's equation in a semiclassical analysis to determine the modified Hamiltonian 
due to the electronic charge reconstruction, but use a time-dependent vector potential to 
describe the electric field that drives current through the system. In the case of cold atoms, 
we discussed nonequilibrium situations corresponding to pulling the lattice in the presence 
of a static trap, rapidly changing the location of the trap and examining the transient re- 
sponse, and changing the interaction strength (say due to a Feshbach resonance, by changing 
the magnetic field) and examining the response to an interaction quench. For the ordered 
phase, we examined the CDW state of the Falicov-Kimball model at zero temperature, with 
an uniform electric field turned on at a specific time. We feel these sum rules will have a 
broad application across a number of different systems and calculations. 

It is obvious that a similar exercise can be carried out for bosonic systems. These will have 
most relevance for cold atom problems, where bosonic atoms are readily available within the 
alkali family. The sum rules for bosons become more complicated, especially at higher orders, 
because we used identities such as ^c^^ = cjq which hold only for fermionic systems. We 
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are currently working on the generalization of these sum rules to bosonic problems which will 
be presented elsewhere. Of particular interest is a Falicov-Kimball model with light fermions 
and heavy bosons, which is realized in K-Rb mixtures used to make dipolar molecules if the 
system is placed on a deep enough optical lattice. 

Acknowledgments 

J. K. F. acknowledges support from the National Science Foundation under grant num- 
ber DMR-0705266. Supercomputer time was provided by the ERDC, ASC, and ARSC 
supercomputer centers of the HPCMP of the DOD. The largest real-axis cold atom calcula- 
tions were performed during a capabilities application project in 2008 on the Cray XT5 at 
ARSC. V.T. acknowledges support by the Department of Energy under grant number DE- 
FG02-07ER15842. We also acknowledge useful conversations with N. Bliimer, M. Kollar, H. 
Krishnamurthy and W. Shen. 



* URL: http://www.physics.georgetown.edu/~jkf 

1 A. B. Harris and R. V. Lange, Phys. Rev. 157, 295 (1967). 

2 J. Hubbard, Proc. Royal Soc. (London) A 276, 238 (1963). 

3 L. M. Falicov and J. C. Kimball, Phys. Rev. Lett. 22, 997 (1969). 

4 W. Nolting, Z. Phys. 255, 25 (1972). 

5 G. Geipel and W. Nolting, Phys. Rev. B 38, 2608 (1988). 

6 W. Nolting and W. Borgiel, Phys. Rev. B 39, 6962 (1989). 

7 M. Potthoff, T. Wegner and W. Nolting, Phys. Rev. B 55, 16132 (1997). 

8 M. Potthoff, T. Herrmann, T. Wegner, and W. Nolting, Phys. Stat. Solidi B 210, 199 (1999). 

9 H. Eskes, A. M. Oles, M. B. J. Meinders, W. Stephan, Phys. Rev. B 50, 17980 (1994). 

P. E. Kornilovitch, Europhys. Lett. 59, 735 (2002). 

1 M. Randeria, R. Sensarma, N. Trivedi, and F.-Ch. Zhang, Phys. Rev. Lett. 95, 137001 (2005). 

2 S. R. White, Phys. Rev. B 44, 4670 (1991). 

3 J. J. Deisz, D. W. Hess, and J. W. Serene, Phys. Rev. B 55, 2089 (1997). 

4 V. M. Turkowski and J. K. Freericks, Phys. Rev. B 73, 075108 (2006); Phys. Rev. B 73 



27 



209902(E) (2006). 

15 V. M. Turkowski and J. K. Freericks, Phys. Rev. B 77, 205102 (2008). 

16 J. K. Freericks, V. M. Turkowski, V. Zlatic, Phys. Rev. Lett. 97 266408 (2006). 

17 V. M. Turkowski and J.K. Freericks, " Nonequilibrium dynamical mean-field theory of strongly 
correlated electrons", in "Strongly Correlated Systems, Coherence and Entanglement", Eds. 
J.M.P Carmelo, J.M.B. Lopes dos Santos, V. Rocha Vieira, and P.D. Sacramento (World Sci- 
entific, Singapore, 2007), pp. 187-210. 

18 J. K. Freericks, Phys. Rev. B 77, 075109 (2008). 

19 J. K. Freericks, V. M. Turkowski, and V. Zlatic, Nonlinear response of strongly correlated 
materials to large electric fields, in Proceedings of the HPCMP Users Group Conference 2006, 
Denver, CO, June 26-29, 2006, edited by D. E. Post (IEEE Computer Society, Los Alamitos, 
CA, 2006), pp. 218-226. 

20 J. E. Hirsch and R. M. Fye, Phys. Rev. Lett. 56, 2521 (1986). 

21 M. Jarrell, Phys. Rev. Lett. 69, 168 (1992). 

22 C. Knecht, N. Bliimer, and P. G. J. van Dongen, Phys. Rev. B 72, 081103(R) (2005) 

23 N. Bliimer, preprint, arXiv: 07 12.1290] (2007). 



24 N. Bliimer, preprint, arXiv:0801.T222\ (2008). 

25 A. N. Rubtsov, V. V. Savkin, and A. I. Lichtenstein, Phys. Rev. B 72, 035122 (2005); P. Werner, 
A. Comanac, L. de Medici, M. Troyer, and A. J. Millis, Phys. Rev. Lett. 97, 076405 (2006). 

26 J. K. Freericks, Transport in multilayered nanostructures: the dynamical mean-field theory ap- 
proach (Imperial College Press, London, 2006). 

27 D. Jaksch, Ch. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, Phys. Rev. Lett. 81, 3108 
(1998). 

28 R. W. Helmes, T. A. Costi, and A. Rosch Phys. Rev. Lett. 101, 066802 (2008) 

29 H. Zenia, J. K. Freericks, H. R. Krishnamurthy, Th. Pruschke, preprint, arXiv:0809.4993 (2008). 

30 S. Okamoto and A. J. Millis, Nature (London) 428, 630 (2004). 

31 J. K. Freericks, V. Zlatic? and A. M. Shvaika, Phys. Rev. B 75, 035133 (2007). 

32 R. E. Peierls, Z. Phys. 80, 763 (1933); A. P. Jauho and J. W. Wilkins, Phys. Rev. B 29, 1919 
(1984). 

33 M. Eckstein and M. Kollar, Phys. Rev. Lett. 100, 120404 (2008). 

34 M. PotthofT, Phys. Rev. B 64, 165114 (2001). 

28 



U. Brandt and C. Mielsch, Z. Phys. B, 75, 365 (1989). 

J. K. Freericks and V. Zlatic, Rev. Mod. Phys. 75, 1333 (2003). 

M. Potthoff and W. Nolting, Phys. Rev. B 59, 2549 (1999). 

Y. Song, R. Wortis, and W. A. Atkinson, Phys. Rev. B 77, 054202 (2008). 

R. W. Helmes, T. A. Costi, and A. Rosch, Phys. Rev. Lett. 100, 056403 (2008). 

O. P. Matveev, A. M. Shvaika, and J. K. Freericks, Phys. Rev. B 77, 035102 (2008). 

W. Shen and J. K. Freericks, unpublished. 

V. M. Turkowski and J. K. Freericks, Phys. Rev. B 71, 085104 (2005). 



29 



